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ABSTRACT 

We estimate the impact of weak lensing by strongly nonlinear cosmological 
structures on the cosmic microwave background. Accurate calculation of large £ 
multipoles requires N-body simulations and ray-tracing schemes with both high 
spatial and temporal resolution. To this end we have developed a new code 
that combines a gravitational Adaptive Particle- Particle, Particle-Mesh (AP3M) 
solver with a weak lensing evaluation routine. The lensing deviations are eval- 
uated while structure evolves during the simulation so that all evolution steps- 
rather than just a few outputs-are used in the lensing computations. The new 
code also includes a ray-tracing procedure that avoids periodicity effects in a 
universe that is modeled as a 3-D torus in the standard way. Results from our 
new simulations are compared with previous ones based on Particle- Mesh simu- 
lations. We also systematically investigate the impact of box volume, resolution, 
and ray-tracing directions on the variance of the computed power spectra. We 
find that a box size of 512h~^ Mpc is sufficient to provide a robust estimate of 
the weak lensing angular power spectrum in the ^-interval (2,000-7,000). For a 
reaslistic cosmological model the power [£{£ + l)C£/27r]^/^ takes on values of a 
few /iK in this interval, which suggests that a future detection is feasible and 
may explain the excess power at high i in the BIMA and CBI observations. 

Subject headings: methods:numerical — cosmic background radiation — cosmol- 
ogy:theory — large-scale structure of the universe 
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Introduction 



Photons from the last scattering surface of the cosmic microwave background (CMB) 
detected at redshift zero, are inevitably lensed by cosmological structures. This lensing 
produces a number of modifications to the CMB angular power spectrum, deviations from 
Gaussianity, and S-polarization. While these effects have been e xtensively studied, and a 



review of their impact can be found in 



Lewis fc Challinorl (120061 ). we are here concerned 



with the lensing due to strongly nonlinear objects such as galaxy clusters and smaller-scale 
structures. In particular, we focus our attention on the estimation of the small angular scale 
(high €) correlations induced by this lensing. This marks a somewhat different direction to 
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N-body simulations are necessary to estimate some CMB anisotropics. This is the 
case for the gravitational anisotropies (weak lensing and Rees-Sciama effects) produced by 
strongly nonlinear cosmological structures. While simulations with sufficient resolution to 
model the formation of galaxies within a n entire Hubble yolume are currently o ut of reach, 



researchers are getting close to this goal (IKim et al. 



2008 



Teyssier et al. 



20091 ). However, 



until such simulations are possible in order to estimate gravitational anisotropies in the 
CMB it is necessary to create an artificial periodic universe where copies of the simulation 
box at various redshifts are replicated back to a given redshift. While numerous alternatives 
exist for how this is achieved all the methods propagate CMB photons along appropriate 
directions within these model universes. When conducting this ray-tracing, box sizes, 
spatial scales and ray directions must be chosen in such a way that the periodicity of the 
artificial universe has no appreciable effects on the resulting anisotropy. The ray-tracing 
procedure used here achieves these goals and was described and applied in earlier papers. 
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where it was also compared with previous methods 1 
Two of these earher papers dealt with weak lensing ( 


Jain et al. 


2000; 


W 


lite & Hu 


2000 


). 


Cerda-Dui 


■an et al. 


2004; 


Anton et al. 



20051 ). and two others with the Rees-Sciama effect (ISaez et al. 



2006 



Puchades et al 



20061) 



The simulations discussed i n these papers were conducted with the Particle-mesh (PM, e.g. 



Hockney fc Eastwood 



19881 ) N-body simulation algorithm. For the large simulation boxes 
required in our study (see below) the spatial resolution of PM simulations is moderate, at 
best. For this reason it was suggested in these papers that better N-body simulations might 
be necessary to estimate good CMB angular power spectra for large multipole indexes 
{£ values). This suggestion is c onfirmed in the p resent study, where AP3M (Adaptive 



Particle-Particle, Particle-Mesh; 



Couchman 



19911 ) simulations with high resolution are 



performed to build up an appropriate periodic universe, in which the CMB photons are 
then moved according to the aforementioned ray-tracing procedure. A CMB angular power 
spectrum has been estimated from i ~ 200 to £ ~ 10000, and the accuracy of this spectrum 
in various ^-intervals is discussed. 

Motivated by a desire to match a number of cosmological constraints, including the 
five-year WMAP (Wilkinson microwave anisotropy probe) observations, high-redshift 
supernovae measurements and galaxy surveys, {e.g. Spergel et a/. 2007; Vianna & Liddle 
1996; Riess et al. 1998; Perlmutter et al. 1999) we consider a cosmological model in which 
the universe is flat and the initial fluctuation spectrum has an inflationary origin. It is 
assumed that only scalar perturbations are present and that weak lensing from gravitational 
waves is negligible. The resulting distribution of perturbations is then Gaussian and the 



perturbations themselves are purely adiabatic. Post inflati on, the power 



1970 



s pectrum o: 



ZeFdovich 



the 



19701 ) 



perturbations is almost exactly of the Harrison-Zel'dovich (iHarrison 
form. This power spectrum is then modified during the radiation-dominated phase, as 
acoustic oscillations in the photon-baryon fluid damp growth, an effect accounted for by 



our use of a CMBFAST calculation of the transfer function. The following cosmological 
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parameters are used for consistency with the standard model (IKomatsu et al.l 12008 



HinshawetaL 



20081 ): (1) a reduced Hubble constant h = 10 '^Hq = 0.7 (where Hq is the 



Hubble constant in units of km s~^ Mpc~^); (2) density parameters = 0.046, Qd = 0.233, 
and = 0.721 for the baryonic, dark matter, and vacuum energy, respectively (the matter 
density parameter is then Qm = + = 0.279); (3) an optical depth r = 0.084 which 
characterizes the reionization; and (4) a parameter cxg = 0.817 normalizing the power 
spectrum of the energy perturbations. 

The layout of this paper is as follows: in Section |5] we outline salient issues in our 
map-making technique. We follow this in Section 12] with an outline of N-body methods 
and our ray-tracing technique. Results from our new simulations are presented in Section 
m along with a critical comparison to earlier work. We then compare our results to current 
observations and conclude with a brief discussion. 



2. Map Construction 

We begin by describing the procedure for constructing lensed maps of the CMB from 
small unlensed maps at the last scattering surface. Note that in this section and those 
following, we choose units such that c = SnG = 1, where c is the speed of light and G the 
gravitational constant. For a quantity A, A^. and Aq denote the value of A at, respectively, 
the time of emission (last scattering surface) and the present. The scale factor is a(t), where 
t is the cosmological time, and its present value, oq, is assumed to be unity, which is always 
possible in fiat universes. The unit vector n defines the observation direction (line of sight). 

Small, unlensed maps of CMB temperature contrasts (A = 6T/T) must be constructed 
to be subsequently deformed by lens ing. These maps have been obtained with a method 



based on the fast Fourier transform (jSaez et al 



19961) and lead to small squared- Gaussian 
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maps of the contrast A. These maps are uniformly pixehsed. The code designed to build 
up the unlensed maps (map making procedure), requires th e CMB angular power spect rum, 



19961) for 



which has been obtained by running the CMBFAST code (ISeljak &: Zaldarriagal 
the model described above. In order to deform the unlensed maps, the lens deviations 
corresponding to a set of directions, covering a n appropriat e region of the sky, must be 



calculated. These deviations are the quantities ( ISeljak 



19961 ): 



Ao 

W{\)V±(j) dX , (1) 

e 

where V_l0 = — n An A V0 is the transverse gradient of the peculiar gravitational potential 
0, and W{X) = (Ae — A)/Ae. The variable A is 



A(a) 



db 



(2) 



In the case of weak lensing, the integral in Eq. ([T]) is usually integrated along straight 
paths, ignoring the small deflections associated with the lensing effect. This approximation, 
commonly known as the Born approximation, is equivalent to keeping first order terms 
in a positional expansion of the transverse potential as a function of the normal ray and 
its lensing of fset. While this approxima tion is known to be comparatively inaccurate at 
small i {e.g. iVan Waerbeke et a/. 1 120011) . at 1000 < ^ < 10000 detailed calculations using 
higher order perturbation theory suggest that corrections to the first or der assumption are 



appro ximately two orders of magnitude lower than the first order term ( IShapiro fc Cooray 



20061). Recent N-bo dy simulation work examining the validity of the Born approximation 



( iHilbert et al. 



20091 ) has shown corrections only begin to become significant at the 5% level 
for I > 20000. Since we are interested in calculating 1000 < i < 10000 values we cautiously 
accept the errors inherent in the first order approach. 



Once the deviations have been calculated, they can be easily used to get the lensed 
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maps from the unlensed ones. This is achieved using the relation 



AM=\in + 6) , (3) 

where and are the temperature contrasts of the lensed and unlensed maps, 
respectively. 

Given the unlensed map A,^, and the map A^ obtained from it after deformation by 
lensing (the lensed map), the chosen power spectrum estimator can be used to get the 
quantities CiiJJ) and Ci{L), whose differences Ci{LU) = Ci{L) — CiiJJ) can be considered 
as an appropriate measure of the weak lensing effect on the CMB. Moreover, a map of 
deformations A^ = A^ — A^, can be obtained and these can be analyzed to get another 
angular power spectrum Ce{D). Since the maps L and U are no t statistically independent, 



Anton et al. 



fl2005h for 



the spectra Ce{D) and Ce{LU) appear to be very different (see 
details). Eight hundred unlensed maps are lensed (deformed) by using the same S field 
and the average Ci{D) and Ce{LU) spectra are calculated and analyzed. Both spectra are 
very distinct measures of the weak lensing effect under consideration. Many times, only 
the customary oscillating Ce{LU) spectrum is shown; however, in a few appropriate cases, 
the Cf(D) spectra are d isplayed. Our power spe ctrum estimator was described in detail in 



Arnau et al. 



(I2OO2I ) and lBurigana fc Saea (120031 ). Results obtai ned with this estirn ator were 



co mpared with those of the code ANAFAST of th e HEA LPix (IGorski et al. 



m 



Arnau et al. 



(120021 ) and also in 



Puchades et al. 



I999I ) package 



(120061 ). These comparisons showed that 



our estimator is a very good one in the case of regularly pixelised squared maps as analyzed 
here. 
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3. N-body simulations and ray-tracing procedure 



3.1. N-body simulation technique 



The primary simulations presented in this work were run using a 



based implementation of the "HYDRA" code (IThacker fc Couchman 



para llel OpenMP- 



20061). This code 



uses the AP3M algorithm to calculate gravitational forces within a simulation containing 
Np particles. In the AP3M algorithm a cubic "base" mesh of size Nc cells per side is 
supplemented by a series of refined-mesh P3M calculations to provid e sub-mesh resolution. 



Grav itational softening is implemented using the S2 softening kernel (IHockney fc Eastwood 



19881 ) which is remarkably similar in sh ape to the cubic spline softening kernel used in many 



treecodes {e.g. 



Hernquist fc Katz 



19891 ). The S2 softening used in the kernel is 2.34 x 5*^ 



where Sp is an equivalent Plummer softening length which we quote throughout the paper 
to enable a simple comparison to other work. The softening length is held constant in 
physical coordinates subject to the resolution not falling below 0.6 of the interpart icle 



Springel et al. 



20051 ) and 



spacing at high redshift. This technique is widely applied {e.g. 
is a compromise between assuring that the potential energy of clusters does not evolve 
significantly at low redshift, while still ensuring structures and linear perturbations at high 
redshift are followed with reasonable accuracy. 

Initial conditions were calcu lated using the standard Zel'dovich approximation 



technique (lEfstathiou et al. 



19851 ). and all simulations were started at a redshift of 2; = 50, 
which is sufficiently early to place modes in the linear regime. To account for the impact 
of varying box sizes, Lbo^., we considered simulations of size 256/i"^ Mpc, 512/i~^ Mpc, and 
1024/i~^ Mpc and a full list of N-body simulation and ray-tracing parameters is given in 
Table [H At z = 6 (the beginning of our ray-tracing epoch in the simulation), the box sizes 
of 256/1-1 Mpc, 512/1-1 y^^^^ 1024/i-i Mpc correspond to square sky patches with 
angular sizes, $map, of 2.48°, 4.96°, and 9.92°, respectively. These are then the sizes of our 
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Table 1. Lensing Simulations 



Name 


Algorithm 




Mode set 


















Direction 












10^" Mq 




/h ^ kpc 






/h~^kpc 






RLS_AA 


AP3M 


512 


A 


512^ 


11 


1024 


12 


512 


6 


25 


0.59 


D512A 


RLS^B 


AP3M 


512 


A 


512^ 


11 


1024 


12 


512 


6 


25 


0.59 


D512B 


RLS 


AP3M 


512 


B 


512-'' 


11 


1024 


12 


512 


6 


25 


0.59 


D612A 


RLS.CA 


AP3M 


512 


C 


512^ 


11 


1024 


12 


512 


6 


25 


0.59 


D512A 


LSJjDA 


AP3M 


1024 


D 


512^ 


88 


1024 


24 


1024 


6 


50 


0.69 


D1024A 


LS31AB 


AP3M 


512 


A 


256^ 


88 


612 


24 


512 


6 


50 


0.69 


D512B 


LS^SIAA 


AP3M 


512 


A 


512^ 


11 


1024 


24 


512 


6 


50 


0.69 


D512A 


LS_AS2AA 


AP3M 


512 


A 


5123 


11 


1024 


36 


512 


6 


75 


0.59 


D512A 


LS^IAA 


AP3M 


512 


A 


512^ 


11 


1024 


12 


512 


6 


12 


0.69 


D512A 


LS^lAB 


AP3M 


512 


A 


512^ 


11 


1024 


12 


512 


6 


12 


0.69 


D512B 


LS^IBA 


AP3M 


512 


B 


512^ 


11 


1024 


12 


512 


6 


12 


0.59 


D512A 


LS^IBB 


AP3M 


512 


B 


512^ 


11 


1024 


12 


512 


6 


12 


0.59 


D512B 


LS^2AA 


AP3M 


512 


A 


512^ 


11 


1024 


12 


512 


6 


40 


0.69 


D612A 


LS_A2AB 


AP3M 


512 


A 


512^ 


11 


1024 


12 


512 


6 


40 


0.69 


D512B 


LS^2BA 


AP3M 


512 


B 


512^ 


11 


1024 


12 


512 


6 


40 


0.59 


D512A 


LS_A2BB 


AP3M 


512 


B 


5123 


11 


1024 


12 


512 


6 


40 


0.59 


D512B 


LSJIBA 


AP3M 


256 


E 


512^ 


1.4 


1024 


6 


512 


6 


15 


0.29 


D256A 


LSJVIGA 


AP3M 


512 


G 


256^ 


88 


612 


24 


512 


6 


60 


0.69 


D612A 


LSJjHA 


AP3M 


512 


H 


128^ 


704 


266 


48 


512 


6 


120 


0.69 


D612A 


PMLS 


PM 


256 


F 


512^ 


1.4 


612 


1000 


256 


6 


500 


0.69 


D266A 



Note. — List of the parameters used in the lensing simulations. Columns three to eight give information associated with the N-body simulations, namely 
the box size, Lj^^^, the set of initial modes (determined by the box size and a single random seed), the number of particles used, Np, the mass of the individual 
particles, Mp, the number of Fourier cells along each side of the simulation, Nc, and the Plummer softening length used, Sp. All N-body simulations were 
started at a rcdshift of z = 50. The following five columns give parameters associated with the lensing calculation, namely the number of rays along the 
edge of the map, N^j^^, the initial rcdshift at which ray tracing begins, z^^, the distance between evaluations on the geodesic Ap^, the angular resolution of 
the CMB lensed maps A.angi and the preferred direction used (see Table 2). 
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constructed CMB maps for each simulation. 



3.2. Ray-tracing technique 



3.2.1. Lensing regimes 



As in our previous work ( lAnton et al 
parts: 



20051 ) we divide the total lensing effect into three 



AWL (A weak lensing), namely the effect due to scales k > 2-n / L^ax (where 
Lmax = 42/i~^ Mpc) at redshifts z <Q. This signal is dominated by strongly nonlinear 
scales 

BWL, the lensing signal due to scales k < 2n/Lmax which corresponds to modes that 
are always in the linear regime down to z = 

CWL, the lensing signal due to scales k > 2ti / L^ax but at redshifs z > Q 



The goal of this paper is to calculate the AWL signal, using ray-tracing through 
N-body simulations, to £ values in the range 1000 — 10000. Since the modes associated with 
the BWL component are specifically chosen to correspond to scales that are linear down 



to z = 0, which is set by wavenumbers k < O.lS/i Mpc ^ (jSmith et al. 



20031), the BWL 



component can be calculated with the linear approach implemented in CMBFAST. The 



CWL component involves modes in the mildly non . 



be evolved via appro ximation schemes (jZel'dovich 



inear regime which can, nonethele ss 



1970 



Sandarin fc Zel'dovich 



1989 



Moutarde et al. 



19911 ). Standard semi- analytical methods designed to study weak lensing 
in the nonlinear regime should also apply in this case, hence, the CWL effect can be 
calculated without resorting to N-body techniques. This calculation can be performed 



using the nonlinear version of CMBFAST which is based on semi-analytical approaches 
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( iLewis fc Challinorll2006l ). For these reasons we begin ray tracing within the simulation 
at z = 6 and consider nonhnear spatial scales having wavenumbers k > 2tt / L^ax- Other 



methods utilizing simulations have also included the contribution to 



above a certain redshift using special techniques (jPas fc Bode 



2008 



lensing from a 


scah 


Carbone et al. 


2008) 



es 



In Figured] we show calculations of the LU angular power spectra for the AWL, BWL 
and CWL lensing regimes for 3200 < ^ < 10000 using the nonlinear lensing implementation 
in CMBFAST. Note that for ^ < 3500 the AWL and CWL effects are too small to be 
calculated with CMBFAST due to the fact that these spectra are obtained as differences 
between other spectra provided by the code. Both the AWL and CWL effects, as calculated 
by CMBFAST, are clearly small with the CWL effect being virtually negligible in the 
3200 < I < 10000 range. As anticipated, the BWL effect is dominant. The calculation 
of the AWL effect provides an initial estimate of the signal we wish to calculate with the 
simulation. We have verified that, for I < 10000, the BWL effects calculated by using both 
the linear and nonlinear methods implemented in CMBFAST are almost identical, which 
proves that the BWL component is linear, as expected. 



3.2.2. Calculation of the lens deviation integral 

An estimate of the lens deviation integral, given in Eq. ([1]), is computed via a numerical 
integration performed along the background null geodesies from the comoving distance 
Dmax — 5900/i~^ Mpc corresponding to 2; = 6 (in the model under consideration) to the 
observer position. The gradient of the peculiar potential V0 used in the integrand of 
Eq. (Il]) is obtained from the simulation but is not exactly the same as that used within 
the N-body calculation. Instead, the gradient is found by subtracting the part of the full 
N-body potential produced by linear spatial scales larger than 42/i~^ Mpc. 
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Fig. 1. — Results from CMBFAST calculations. Top-left, top-right and bottom left panels 
exhibits the LU angular power spectrum for i > 3200 in the cases AWL, BWL, and CWL, 
respectively. Bottom-right panel is the LU spectrum corresponding to the BWL effect for 
£ < 10000. This effect dominates against AWL and CWL for £ < 3000. 



Specifically, our algorithm for determining the potential gradient, which is proportional 
to the force, is as follows: 
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1. Decide upon the direction of the normal rays representing the geodesies (see 
Section 

2. Assuming the Born approximation and using the photon step distance Apg (see 
Section I3.2.3P determine all the evaluation positions and times on the geodesies within 
the simulation volume from z = 6 down to the final redshift 

3. Associate test particles with each of these positions and times 

4. At each time-step of the N-body simulation (while it is running) determine which test 
particles require force evaluations 

5. At each test particle position evaluate the force on the test particle using the 
long-range FFT component and short-range PP correction as in the HYDRA 
algorithm 

6. During the FFT convolution for the test particles eliminate contributions from scales 
larger than A2h^^ Mpc by removing the signal from wavenumbers satisfying k < 0.15h 
Mpc"^ 

7. If the evaluation time for a point on the geodesic lies between two time-steps calculate 
a linear interpolation of the two forces from the time-steps that straddle the correct 
time 

8. Resolve the force into its transverse component and hence recover the transverse 
component of the potential gradient 

Once all the potential gradients are evaluated we can calculate the lensing deviation 
integrals. We emphasize that our algorithm ensures the potential gradient along normal rays 
is calculated very accurately and corresponds exactly (modulo the removal of power from 
scales greater than 42h~^ Mpc) to that in the simulation. We do not resort to smoothing 
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onto grids or the creation of lensing planes "on the fly" . The main drawback of our method 
is that since lensing is now an integral part of the N-body simulation, a different preferred 
direction requires a new N-body simulation. 



3.2.3. Photon propagation paths 



To calculate the photon propagation, the CMB photons are moved through the 
simulation volume along specially chosen paths to avoid repeatedly sampling the same 
structures. This approach uses all steps from the simulation and by using the periodicity 
of the box volume there are no discontinuities in the matter field anywh ere. In princi ple 



this a pproach is simi 



2000 



Hamana et al. 



ar to tiling methods t hat are used elsewhere {e.g. 



2001 



Sato et al. 



White fc Hu 



20091 ). These methods are based on independent 



PM simulations with decreasing sizes that telescope in resolution along the line of sight. 
In our "tiling" the matter field is constantly being updated by the AP3M code as the 
time-step changes and, consequently, it is not limited by the box size. We also take care 
to ensure paths are taken which avoid, as much as possible, periodicity effects. This is 
notably different fro m other approaches using random translations and orientations {e.g. 



Carbone et al. 



20081 ) of the simulation volume which, unavoidably, have discontinuities 



at adjoining radial shells. While the signal from such discontinuities is likely small, our 
method has the advantage of avoiding it completely. 

Choosing the directions for the ray propagation is not entirely trivial. For the directions 
parallel to the box edges, periodicity effects are clearly very strong. Photons moving 
along these directions would pass close to the same structures in successive boxes and, 
consequently, the lensing effect of the same structure would be included a number of times 
(one per crossed box). This repetitio n would lead t o a fa lse magnification of len s deviations. 



However, as has been emphasized in 



Anton et al 



fl2005[ ) and 



Saez et al. 



(|2006[ ). periodicity 
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effects are negligible along certain directions, hereafter called preferred directions. 

In order to define preferred directions, it can be assumed that: (i) the x, y, and z 
axes are parallel to the box edges; (ii) the angles 9 and ip are spherical coordinates defined 
with respect to these axes; (iii) photons moving along the direction n, cross the {y, z) face 
of a box at point P, and the next (y, z) face at point Q; and (iv) if the segment PQ is 
projected onto the {y,z) plane, the length of the resulting projection is Cpq- Then, if the 
condition C^p^ > Lmax is satisfied, the direction n is assumed to be a preferred one. Taking 
into account this definition of preferred directions and the cutoff performed at the scale 
Lmaxi it is clear that photons moving along preferred directions enter successive boxes 
through independent uncorrelated regions. Moreover, for these directions and the box sizes 
of interest, it can be easily verified that the CMB photons can travel from z = 6 to z = 
(~ 5900/i^^ Mpc) through different uncorrelated regions (without repetitions) located in 
successive simulation boxes. Periodicity effects can, therefore, be assumed to be negligible. 
The total number of crossed boxes is denoted by N^r- 

Some preferred directions used in this paper are defined in Table [2], where the first 
column gives the names of all these directions. Each of them has been used in universes 
covered by simulation boxes with a given size, as detailed in the second column. The third 
and fourth columns show, respectively, the angles 6 and if (spherical coordinates) defining 
the corresponding direction. The distance Cpq defined in the previous paragraph is given in 
the fifth column and, finally, in the last column the number of boxes crossed by the CMB 
photons from redshift z = 6 to the present is given. From the values of Table [21 it is clear 
that the distance (p^ is much greater than Lmax = 42/i~^ Mpc in all four cases, which 
ensures that other directions close enough to the listed ones are also preferred directions. 
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Table 2. Preferred directions for ray tracing. 



Direction 


Li box 


e 


95 






D256A 


256 


76.76 


11.31 


79.98 


22.00 


D512A 


512 


76.76 


11.31 


159.95 


11.00 


D512B 


512 


68.43 


18.435 


273.20 


10.17 


D1024A 


1024 


59.19 


26.57 


853.33 


4.43 



Note. — Distances Li,ox and are given in units 
of Mpc, and the spherical coordinates 9 and (/? 
in degrees 
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We proceed as follows: any direction of Table [2] is assumed to point toward the center 
of a squared map, whose angular size corresponds to the box size appearing in the second 
column (as noted above). Then, these squared maps are uniformly pixelised by choosing a 
certain number of pixels, Npi^, per edge. The angular resolution is then /S.ang = ^map/Npix. 
The directions of all the pixels are preferred ones and, consequently, lens deviations can be 
calculated for each pixel — with no significant periodic effects — across the full map. 

The parameters involved in the ray-tracing procedure are thus summarized: a number 
of directions, N^ir, per edge of the squared CMB map (one per pixel, Ndir = Npi^); sua initial 
redshift, z^, for the calculation of lens deviations; a step, Ap^, to perform the integral in 
Eq-dl]) (hereafter called the photon step); and the angles 9 and ip defining the preferred 
direction. 

Hereafter, a lensing simulation (LS) is the calculation of the (5-deviations along the 
pixel directions, plus the construction of the A^, A^, A^ maps and the estimation of 
Ci{D) and Ce{LU) (angular power spectra). An LS is characterized by the parameters 
and initial conditions required by the N-body simulation together with the parameters of 
the ray-tracing procedure. The lensing simulations (LSs) obtained from the parameters: 
Lbox = 512/1-1 Mpc, Np = 512^, = 1024, Sp = I2h~^ kpc, N^^r = 512, Zir^ = 6 (as in 
all the simulations used in the paper), and Aps = 25h~^ kpc, are hereafter called reference 
lensing simulations (RLSs). The angular resolution of these simulations is Aang — 0.59' 
[i ^ 18,600). There are an infinite number of possible realizations of this type of LSs 
corresponding to different initial conditions for the N-body simulation as well as to distinct 
preferred directions. We also consider the effect on the LSs of changing parameters to ensure 
that the calculation of the power spectra Ce{D) and Ci{LU) is as robust and accurate as 
possible. 
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4. Results 

4.1. Effect of changing the final ray-tracing redshift 

To determine the impact of changing the final ray-tracing epoch, the RLS_AA lensing 
simulation (based on the preferred direction D512A, see tables |2] &[!]), was used to estimate 
lensing from = 6 to final redshifts between 0.5 and 0. Results are shown in Fig. [21 where 
the D angular power spectra corresponding to the final redshifts 0.3 (dotted line), 0.2 (solid 
line), and 0.0 (dashed line) are presented. Clearly, these three spectra are very similar. We 
find the same result for the LU spectra (these results are omitted since the addition of five 
oscillating lines would lead to a confusing plot and add little information). Based upon 
these results we conclude that the signal produced between redshifts 0.2 and contributes 
negligibly to the total lensing. Given this fact, and that the CPU time required to evolve 
between z = 0.2 and 2; = is comparatively lengthy for our code (due to the absence of 
individual particle time-steps), our calculations are performed between redshifts Zj„ = 6 
and Zend = 0.2. The same study has been done for other RLSs corresponding to different 
initial conditions and preferred directions. The conclusions are the same in all the cases, 
namely a negligible effect be tween z = 0.2 and z = 0. We note that this conclusion is in 



close agreement with that of 



Carbone et al. 



(I2OO8I ) who showed that stopping at z = 0.22 



produced a deficit in the low £ signal but a negligible difference for £ > 350. 



4.2. AP3M simulations contrasted to PM 

We next compare LSs constructed from PM and AP3M N-body codes to examine the 
impact of sub-mesh-scale resolution on the lensing signal. In Fig. [3] we plot the Ci{LU) 
(top panel) and Cf {D) (bottom panel) for RLS_BA (solid line) along with the same spectra 
for PMLS (dashed line), obtained from a PM code. Details of the PM code can be found 
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Fig. 2. — D angular power spectra corresponding to RLS_AA at different final redshifts. 
Curves corresponding to z^nd values of 0.5 (triple-dot-dash), 0.4 (dot-dash), 0.3 (dots), 0.2 
(solid), and 0.0 (dashes) are shown. 
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in iQuilis et al.l (Il998l ): this code has been used in previous papers {e.g. lAnton et aLll2005l ). 
We emphasize that the same ray-tracing procedure, described in Section [31 is used in both 
codes. The parameters of the PMLS (see Table [T]) are as follows: L^ox = 256h~^ Mpc, 
Np = 512^, Nc = 512 , Ndir = 256, and Aps = 0.5h~^ Mpc, and the preferred direction is 
D256A (see Table [2]). The effective resolution Eres of the PM code is two cells (1 Mpc 
in this case); hence, in the PMLS, the photon step Apg must be smaller than Ih^^ Mpc to 
take advantage of the PMLS resolution. It has been verified that Ap^ = 0.5h~^ Mpc (half 
of Eres) is a good value for the photon step (smaller valu es lead to very similar results). In 



keeping with other N-body work {e.g. Moore et a/. 1 119981 ) , the true effective resolution of 



the AP3M simulations is estimated to be Eres ~ 55*^, which means that, in our RLSs, the 
photon step is a little smaller than a half of the effective resolution. 

Examination of Fig. 12] shows that the PMLS traces the peaks, but the amplitudes 
are too small and, moreover, the Ce{LU) quantities quickly tend to zero as i increases. 
However, in the AP3M case, the peaks have greater amplitudes and, furthermore, a signal 
of a few micro-Kelvin appears for high i values. Similar conclusio ns follow from the D 



Anton et al. 



(120051) 



spectra of the bottom panel, which can be directly compared with 
where this kind of spectrum was used. These results show that the PMLS underestimates 
the lensing signal we are calculating. Hence, high resolution N-body simulations, be they 
AP3M or computed using some alternative algorithm, are necessary to estimate CMB weak 
lensing from strongly nonlinear structures. 



4.3. Variance of power spectra due to modes in the initial conditions 

Having shown that the results from the AP3M simulations are distinctly different from 
those of the PM simulations, we must still account for variability of results due to changing 
the random modes in the initial conditions of the simulation. To this end we show, in 
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Fig. 3. — Top: LU angular power spectra corresponding to the RLS_BA (solid line) LS, 
and to the PMLS LS defined in the text (dashed line). Bottom: D spectra for the same 
simulations as in top panel. 
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the top panel of Fig. HI the LU angular power spectra corresponding to three distinct 
RLSs, namely, RLS_AA, RLS_BA & RLS_CA, each drawn from a different realization of 
the modes but using the same preferred direction {Dbl2A in Table [2]). The results are 
clearly qualitatively very similar and the only noticeable quantitative differences are below 
I ~ 1000. This is not entirely surprisin g since the lower the I value the larger the sample 



variance (IKnox 



1995 



Scott et al.l 119941) ■ Additionally, the largest deviation occurs at the 



point where the signal is weakest, again an unsurprising result. 

These results suggest that a good average LU spectrum can be achieved from just 
a few chosen RLSs simulations, perhaps as few as two. We have verified this assertion in 
the bottom panel of Fig. HJ by plotting the average of two of the RLSs as compared to the 
average of all three RLSs. Both average spectra are so similar that only two RLSs suffice to 
get a very good average LU spectrum in the ^-interval under consideration. Unsurprisingly, 
the average signals for £ < 1000 are in good agreement, as would be expected if the 
differences in the single realizations were due to sample variance. For ^ > 1000, the spectra 
of the different realizations are so similar that a unique RLS leads to a rather good LU 
spectrum of the lensing effect under consideration (L < L^ax and z < Zin = 6). 

We emphasize that the sample variance we have been talking about is inherent in the 
lensing calculation itself, rather than the CMB maps themselves. For each RLS, the spectra 
are calculated from 800 small maps having a size of 4.96 degrees; which is a necessary step 
to understand small sample variances (small deviations from spectrum to spectrum). 



4.4. Variance of power spectra due to different preferred directions 

In terms of the variance produced by different preferred directions, the ergodic nature 
of the ray-tracing means that we might reasonably expect the same level of variance seen 
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Fig. 4. — Top: Three LU angular power spectra extracted from distinct RLSs, namely 
RLS_AA (solid), RLS_BA (dashed), RLS_CA (dotted). Bottom: the dotted hne gives the 
average LU spectrum of the three RLSs of the top panel, whereas the solid line corresponds 
to the average of two spectra, namely just RLS_BA and RLS_CA. 
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for the different N-body realizations with differing modes. Nevertheless, choosing a different 
preferred direction for a given RLS means that the photons move through different regions 
of the simulation box and cross box faces in different places and a different number of 
times. So the trajectories are not fully equivalent, in turn producing lensing deviations that 
are subtly distinct. To evaluate this difference we plot, in Fig. [5l two LU angular power 
spectra for the RLSs RLS_AA and RLS_AB which use the same initial modes, but have 
preferred directions D512A and D512B. The similarity of the two curves indicates that, as 
anticipated, results are almost independent of the chosen preferred direction, and that our 
ray-tracing procedure is robust. 

4.5. Impact of simulation box size on the power spectra 

Thus far we have considered AP3M simulations with box sizes Lhox = 512h~^ Mpc, 
and hence the LU and D spectra have been obtained from small maps of angular size 
^map = 4.96°. We next examine whether these values are too small to make robust 
predictions. 

In order to answer this question, we conducted an LS, labeled LS_LDA in Table [H with 
Lbox = I024:h~^ Mpc and N^ir = 1024, which leads to a lensed map of size $map = 9.92° and 
with the same angular resolution as in the RLSs. Since simulating a larger volume at fixed 
particle resolution requires a larger softening and the corresponding photon step, the values 
Sp = 24/i"^ kpc, and Apg = 60h^^ kpc were used. The preferred direction was D1024A 
(see Table [2]). To compare to this simulation we also ran another LS, denoted LS_MAB, 
with parameters set to mimic the LS_LDA simulation but in l/8th the volume and using 
l/8th the number of particles. The LU spectra obtained for the two simulations are shown 
in Fig. [6l and are extremely similar. This is despite the fact that these two simulations 
have differing initial conditions, different box sizes, and different numbers of crossed boxes. 
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Fig. 5. — LU angular power spectra obtained from RLS_AA and RLS_AB which share the 
same N-body initial conditions but have distinct preferred directions. The dashed (solid) 
line corresponds to RLS_AA (RLS_AB) and hence the direction D512A (D512B). 
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Ncr- Hence we conclude that a box size of 512/i~^ Mpc — the same as that of the RLSs — is 
large enough to get very good angular power spectra for the range of ^ considered. Even 
smaller sizes, e.g. , 256/i"^ Mpc, can be used when only large enough £ values are under 
consideration (see below). 



4.6. Noise in the power spectra 



Although we do not directly calculate the 2d lensin g potential or converg ence we can 



Carbone et al. 



mm . At our 



examine the power spectrum of the angle a = as in 
final redshift of 0.2, the angular resolution of our RLSs (0.59') corresponds to a comoving 
separation of lOO/i"^ kpc, which is still larger than our nominal N-body resolution of 
bSp = 60h^^ kpc. Since this is the effective resolution in the ray-tracing we do not need to 
worry about the intrinsic resolution of our maps falling below that of the N-body simulation. 

In Fig. [7] we plot the power spectrum of the lens deviations in the interval 
2000 < £ < 10000. We have quantified the noise in this signal by considering box sizes 
of 1024/i~^ Mpc, 512h~^ Mpc, 256h~^ Mpc although our smallest box has an effective 
angular resolution that is twice as small as the other two simulations. For comparison 
to semi-analytic methods we also include the total lensing signal predicted by CAMB in 
this £ range. Comparing to a running average constructed by binning the average signal 
in £ ± 100, the RMS deviations relative to the overall signal are 4.1%, 5.3% and 6.2% for 
the three box sizes, respectively. Although we do see an increase in the relative noise with 
reducing box size, the overall noise is much less than the underlying signal. Further, the 
absolute value of the signals are in close agreement as well, the largest difference being 
~ 15% at £ = 2000 and the qualitative shape of the power spectra agree well at least for 
i < 7000. For i = 10000 we see a difference of 50% in the signals, and we are not confident 
of convergence in the interval 7000 < £ < 10000 (see Sections 14. 8[ 14.91 for an extensive 
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Fig. 6. — LU angular power spectra extracted from two LS simulations (LS_LDA, solid 
line and LS_MAB, dashed line) that differ in box sizes but use similar effective particle 
and map resolutions. The simulations share the same softening and photon steps, namely 
Sp = 2Ah-^ kpc and Ap, = QOh-^ kpc, but for the LS_LDA LS Uox = 1024/1"^ Mpc, 
Np = 512^ while for the LS_MAB LS, Lbox = 512h-^ Mpc, Np = 256^. The initial conditions 
of these LSs are fully independent and different preferred directions are used. 



- 28 - 



discussion of convergence issues). 

4.7. Impact of spatial resolution 

Having examined sampling issues related to box sizes and preferred directions we 
now address the impact of force softening and the step size in the ray-tracing. As a first 
investigation, we varied parameters Sp and Aps, in lock-step with the remaining RLS 
parameters held fixed. A ratio Eres/^ps — 2.5, identical to that of the RLS simulations 
has been assumed in all the cases, thus maintaining a constant ratio between the photon 
step and the effective resolution {i.e. 2.5 photon steps per effective resolution interval). We 
considered the following new pairs of Sp and Ap^: (i) Sp = 24:h~^ kpc and Aps = 50h~^ kpc 
and, (ii) Sp = 36h^^ kpc and Aps = 75h~^ kpc, and call the LSs associated with these 
parameters LS_AS1AA and LS_AS2AA respectively. The C^{LU) for these LSs, along 
with RLS_BA for comparison, are given in Fig. [HI Up to £ ~ 7000 the three spectra are 
very similar, although for I > 7000 we see a small separation in the results. For i > 5000, 
the trend with decreasing softening and photon-step appears to be systematic, a smaller 
softening leading to smaller Ci{LU) coefficients. Nevertheless, a smaller softening does not 
necessarily mean a more realistic simulation and, consequently, we cannot decide which of 
the three lines of Fig. [8] is closer to the true spectrum. Regardless of this issue, any of these 
lines gives a good estimate of the required spectrum for £ < 7000. 

To address the role of the photon step alone, we next held fixed all the RLS parameters 
except the photon step, Ap^. Along with our original choice of A^^ = 25h~^ kpc, two new 
values were considered: Ap^ = 12h~^ kpc (~ Eres/5, corresponding to LSs LS_A1AA, 
LS_A1AB, LS_A1BA, LS_A1BB), and Ap, = 40/i"^ kpc (a little smaller than ~ Eres, 
corresponding to LSs LS_A2AA, LS_A2AB, LS_A2BA, LS_A2BB). For each of the Ap^ two 
distinct preferred directions, namely D512A and D512B, were considered along with two 



-29- 




oE , , , I , , , I , , , I , , , : 

2000 4000 6000 8000 10000 

I 

Fig. 7. — Top, medium and bottom panels correspond to box sizes of 1024/i~^ Mpc 
(LS_LDA), 512/1-1 Mpc (RLS_BA ), and 256/1-^ Mpc (LS_HEA). Dotted lines correspond 
to the correlations extracted from the simulated a maps. Solid lines are averaged spectra 
in running bins of size i. ± 100. The mean amplitude and the rms of the relative errors 
(dotted-solid/solid) are (-0.00072,0.041), (-0.00065,0.053), and (-0.00056,0.062) in the top, 
medium, and bottom panels, respectively. The dashed line in the bottom panel corresponds 
to the total lensing signal predicted by CAMB. Note that the solid and dotted lines in all 
panels correspond to the AWL lensing component. 
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Fig. 8. — LU angular power spectra extracted for RLS_BA (solid line), and from other two 
LSs with the same parameters and initial conditions, but differing softenings and photon steps 

where, Apg = 50h^^ kpc and Sp = 24/i^^ kpc (LS_AS1AA, dotted line), and A^^ = 75h~^ kpc 
and Sp = kpc (LS_AS2AA, dashed hne). 
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distinct sets of initial modes, giving rise to four LSs. The four spectra for each Aps value 
were then averaged and the resulting spectra plotted in Fig. O Note, that since the error 
in the lensing deviation integrals is reduced with decreasing Apg, we can reasonably argue 
that our best estimate of the lensing signal is given by the minimum Aps- We see that in all 
three cases the Ce{LU) quantities are in close agreement up to ^ ~ 7000. In the ^-interval 
[7,000-10,000], the separations between the three curves increase, with a systematic trend of 
a lower signal with decreasing photon-step. Since all these spectra have the same simulation 
softening it is reasonable to conclude that the trend we observed when comparing different 
Sp and Aps may actually be a function of the change in Apg. Given that the smallest Aps 
exhibits the lowest signal, we thus conclude that our best estimate of the signal actually 
has the smallest correlations for i > 7000 values. Nonetheless, for i < 7000, the RLS closely 
follows our best estimate, hence for i < 7000 we are very confident that the RLS leads 
to a robust estimate of the lensing effect we are studying. Even if the Aps = 12h^^ kpc 
value leads to the most accurate spectrum (which is not certain), the computational cost of 
constructing these LSs is significantly greater than that of the RLSs, while only providing a 
mild improvement in accuracy for i > 7000. 

We have also attempted to reproduce results similar to that found elsewhere using 
interpolation methods. We ran another simulation using RLS_AA parameters, except that 
this time the deflections were calculated using the 8 nearest neighbours. This keeps the 
mass resolution of the simulation fixed and does not alter the gravitational calculation, 
while at the same time lowering the effective resolution of the lensing methodology. The 
resulting power spectrum is shown in Fig. ITOland shows a decay ing signal at high i which 



is similar to that found in earlier work {e.g. 



Das fc Bodd (120081 )). This shows that as we 



degrade the resolution of ray-tracing method we do indeed recover previous results. 
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Fig. 9. — Average LU angular power spectra extracted from RLSs (solid line), and from other 
two LSs with the same parameters, excepting the photon step Ap^, whose values are Apg = 
40h~^ kpc (dotted line) and Apg — 12h~^ kpc (dashed line). The dashed, solid, and dotted 
lines corresponds to the smallest, medium (RLS), and greatest Aps values, respectively. 
Averages have been performed by using four appropriate simulations (see text). 
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Fig. 10. — LU angular power spectra for the RLS_AA simulation (solid line) as compared 
to the same simulation but where deflections are calculated by including an average over 
the 8 nearest geodesies (dotted hne). This reduces the resolution of geodesic method, but 
maintains the same resolution in the gravitational solver. Also shown (dashed line) is the 
highest resolution (but small volume) LS_HEA simulation. As shown in more detail in 
figure 11, higher resolution reduces the impact of discreteness effects, which in turn reduces 
the asymptotic high ^ signal. However, sample variance issues in the LS_HEA simulation 
preclude us for drawing firm conclusions at present. 
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4.8. Impact of mass resolution on convergence 

Since our ray-tracing method samples the local gravitational field on scales down to a 
few tens of kpc, it is important to investigate the role of mass resolution on our simulation. 
While we have investigated box size earlier, which at a fixed resolution leads to a change in 
mass resolution, in this section we focus on changing the mass resolution in a fixed volume, 
and also include one additional simulation in a smaller volume at our maximum mass 
resolution. 

We ran additional simulations with 128^ and 256^ particles in the RLS box size, 
namely 512h~^ Mpc. These simulations are labeled LS_LHA and LS_MGA in Table [TJ 
Notable other parameters were Sp = A8h'^^{2Ah^^) kpc, Ap^ = 120h^^{Q0h^^) kpc and 
particle masses were 7.0 X 10^2 (^8 3 X iQii) Mq respectively. In Fig. [H] we plot the power 
spectrum results for these simulations, along with an RLS, another 512^ simulation with 
a shorter photon step (LS_1AB), and one final simulation labeled LS_HEA that uses a 
volume l/8th the size of an RLS but with the same number of particles. The parameters 
of this LS were Lhox = 256h~^ Mpc, Sp = 6h^^ kpc and Ap^ = 15h^^ kpc. Note the two 
512^ simulations have the same mass resolution but they bracket a range of photon-step 
values, a full estimate of the uncertainty in the signal at this mass resolution is given in 
Figure 13, where the spectra of these two cases appear inside the band of uncertainty. Due 
to the different box size of the LS_HEA run we could not maintain the same modes across 
all runs, and hence we used different initial modes at different resolutions. Taken together, 
these simulations span a ratio in mass resolution of 512. 

The plot shows that when PP corrections are included there is a resolution-dependent 
up-turn in the power spectrum for i > 5000. Impr oving the mass r esolution reduces the 



impact of N-body discreteness on the lensing {e.g. I Jain et al.l 120001 ) to the point where 



the highest resolution simulation (LS_HEA) appears to be free of this contamination for 
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Fig. 11. — Comparison of convergence in the LU angular power spectra for different mass 
resolutions. The dot-dashed line corresponds to the 128^ simulation LS_LHA, the dashed 
line to the 256^ simulation LS_MGA, the triple-dot-dashed line to RLS_AA, the solid line 
to LS_A1AB, and the dotted line to LS_HEA. The simulation boxes are 512/i~^ Mpc except 
for the LS_HEA simulation which uses l/8th the volume. The two 512'^ simulations have 
the same mass resolution but they bracket a range of photon-step values (see Table [T] for 
simulation details), a full estimate of the uncertainty in the signal at this mass resolution is 
given in Figure 13. 
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Fig. 12. — LU angular power spectra extracted from two LS obtained from fully independent 
initial conditions, but differing in short-scale resolution. The dotted line corresponds to 
LS_HEA with parameters L,,ox = 256/1"^ Mpc, Np = 512^, A^^ = 1024, N^ir = 512, Sp = 
6 X 10~^h~^ Mpc, and A^^ = lbh~^ kpc. The sohd line corresponds to LS_A1AB which has 
Aps — 12 kpc and the remaining parameters identical to those of the RLSs. 
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£ < lO'^. Turning to lower I values, the peak of the power spectrum around £ ~ 2000 
appears to be converged, modulo the variation due to the initial modes, for resolutions of 
512^ and greater. The results for ^ < 1000 also appear to show errors consistent with the 
variation in initial modes observed in Section I4.3[ indicating results are converged in this 
range for all resolutions. 

4.9. Is convergence at £ > 7000 possible? 

To address the question of the correct form of the spectrum at £ > 7000 we utilize 
the RLSs and the LS_HEA simulation. The effective resolution Ey.f.s — 30/i^^ kpc of the 
LS_HEA simulation is half of that attained in the RLSs and thus small scales should be 
represented more accurately. As a result of the reduced box size, the angular size of the 
lensed maps is $map = 2.48° and the angular resolution limit ^ang = 0.29' (£ ~ 37000) is 
half that of the RLSs {i.e. higher resolution). 

The resulting spectrum is plotted in Fig. [12] as a dotted line, and compared to our 
previous best estimate for the signal at £ > 7000, namely that from LS_A1AB that was 
constructed with Ap^ = 12h^^ kpc and the remaining parameters identical to those of the 
RLSs (the preferred direction is D512B). As expected, given that the new spectrum is 
obtained from smaller 2.48° x 2.48° maps, some uncertainties in the C^{LU) quantities are 
observed for £ < 2000. Note that for the 4.96° x 4.96° maps used in the RLS we expect 
uncertainties for i < 1000. However, for C. > 2000, the spectrum of the smaller box (dotted 
line) is expected to be more accurate than our previous LS. Although we again observe 
a systematic trend of better resolution producing a lower signal at high i > 7000, the 
results are nonetheless in good quantitative agreement for 2000 < i < 7000. We can hence 
be confident that our RLSs give also a rough but useful estimate of the angular power 
spectrum for 2000 < i < 7000 (supporting our assertion in the previous section). Finally, 
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for 7000 < £ < 10, 000, the smaller box produces a decrease, whereas the previous LS 
increases slightly. This unfortunately suggests that the N-body simulations and the CMB 
lensed maps used to derive the RLSs are not good enough for £ > 7000. We hence conclude 
that future simulations with higher resolutions are clearly necessary to go beyond £ = 7000. 



5. Discussion and observational implications 

To bring together all that we have learned thus far, in Fig. [T2] we plot all the LU 
angular power spectra from our simulations (excepting simulation LS_LHA with Np = 128^ 
particles) in the log(£)-interval (3.5, 4.0), which we reemphasize corresponds to the AWL 
effect described earlier (scales smaller than 42h~^ Mpc at z < 6). The AWL and the 
complementary lensing effect BWL estimated with CMBFAST are also included in the 
same plot. Since we demonstrated the CWL component is negligible in the interval under 
consideration (see the bottom left panel of Figd]), it is not represented. The lower solid line 
displays the same spectrum as the dotted line of Fig. [121 which corresponds to the LS_HEA 
LS, which has the smallest value of Sp and the highest mass resolution considered (eight 
times better than that of the RLSs). The upper solid line shows the same spectrum as the 
dashed line of Fig. [HI which corresponds to the LS_AS2AA LS, which has the same mass 
resolution as the RLSs but an Sp value 3 times larger (and 6 times larger than that used 
in the LS_HEA LS). Despite the large differences in the numerical parameters defining the 
N-body simulations, we see that from £ ~ 3000 to ^ ^ 7000, all the spectra lie in a narrow 
band (bounded by the solid lines of Fig. [T3|) . having a width close to half a micro-Kelvin. 
In addition to the simulation lines (solid and dotted) we have plotted estimates of the 
l ensing signal based upon semi-analytical approaches implemented in the CMBFAST code 



( [Lewis fc Challinor 



20061 ). By using the nonlinear version of this code with a minimum scale 



of 30/i ^ kpc (as in our LS, leading to the lower solid line of Fig. [13]), we have estimated 
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the AWL (dash-dots hne in Fig. [13] and top left panel of Fig. [T]) and the BWL components 
(dashed line in Fig. [T3] and top right panel of Fig. [T]) of the weak lensing. Finally, the 
CMBFAST code has been also used to calculate the CMB angular spectrum, in the absence 
of lensing, for the model described in Section |2] (dash-three-dots line). 

The complementary effect BWL, not calculated within the simulation, is due to 
scales greater than Lmax = 42/i~^ Mpc. While the absence of this power has significant 
consequences for small i values we will now show (with CMBFAST) that in the interval we 
are examining it is small. If the scale Lmax = 42/i^^ Mpc is placed on the last scattering 
surface at [z ^ 1100), it subtends an angle of 14'. 63. However, in Fig. [T3]we are plotting 
i values corresponding to angles of ~ 3'. 4 and ~ 1', i.e. considerably smaller than the 
scale set by Lmax on the last scattering surface. By this reason, the power of the dashed 
line, which is due to the linear scales greater than Lmax, is expected to be small in the 
log(£)-interval of Fig. [131 CMBFAST estimates have confirmed this smallness (dashed line 
of Fig. [T3|) . Nevertheless, the BWL effect is very important for much smaller i values (see 
bottom right panel of Fig. [1]). 

The complementary CWL effect (see the bottom left panel of Fig. [1]) is too small to 
be significant in the log(£)-interval (3.5, 4.0). The AWL effect obtained with CMBFAST 
(nonlinear version) is shown in the dot-dashed line of Fig. [131 This effect appears to be 
smaller than the AWL effect numerically estimated here (with AP3M simulations), whose 
power is found to be close to 2 for 3000 < £ < 7000. The total lensing would be 
obtained by adding this power and that of the dashed line {BWL component). We thus see 
that the total lensing effect has a power of a few /i K. 



b compare this result with other simulation work, we turn to the work by 



Das fc Bode 



(I2OO8I ). Since our map-making method does not allow us to plot the lensed and 



unlensed Ci we have converted their results (specifically those in their Figure 4) to the 
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Fig. 13. — Angular power spectra, in fiK, as functions of log(£). Solid and dotted lines 
represent all the spectra displayed in Figs. HlfT2l The lower solid line displays the same 
spectrum as the dotted line of Fig. [121 which corresponds to the LS_HEA LS, which is the 
highest resolution simulation considered. The upper solid line shows the same spectrum as 
the dashed line of Fig. El which corresponds to the LS_AS2AA LS, which has the same mass 
resolution as the RLSs but an Sp value 3 times larger. The long dashed line corresponds 
to the geodesically averaged signal (mimicing reduced resolution), as shown in Figure 10. 
The dash-three-dot line is the CMB spectrum in the absence of lensing, and the dash-dot 
(dashes) line is the AWL (BWL) effect obtained with the CMBFAST code in the case of 
nonlinear lensing including structures with sizes L > 30h~^ Kpc. 
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+ l)Ci{LU) / {27f)Y^'^ values we use in our paper . In performing this conversion we have 
multiphed their power spectra by a factor of 2.90 x 10^'^ fiK^ which includes the effective 
fractional sky coverage, We hence find values of 1.1 jjK at £ = 4000, 0.7 /iK at 

i = 5000, 0.6 /iK aX i = 6000 and 0.4 /iK at i = 7000. This compares with our own signal 
of 2.0 ± 0.4 /iK across the entire range. However, we caution against over-interpreting 
this comparison. Our cosmolo gy and assurnptions about reionization and backgrounds are 



slightly different from those in 



Das &: Bodd ( 120081 ). however it is clear that we find a higher 



signal at ^ > 1000 when short range corrections are included. 

As it is apparent from Fig. [121 the AWL effect clearly dominates the primary 
anisotropics for i > 4200 {i.e. the continuous and dotted lines are well above the dash-three- 
dots one). However, the Sunyaev-Zel'dovich (SZ) effect is also important at these scales, and 
foregrounds must be pr operly accounted for . The Berkeley-Illinois-Maryland Association 
(BIMA), has reported (jPawson et a/.ll2006l ) a temperature power + l)C*^/27r]^/^) of 



14.881^88 f^K at 68 % con fidence, 



one must have as — 1.03 fiCooray 



or i = 5237. In order to ex plain 14.88 fiK as an SZ effect, 



2002 



Dawson et al. 



20061 ). but this is in comparatively 



strong disagreement with data from WMAP (see Section [T]) that l ead to as — 0.82. Hence 



taking into account that the SZ temperature power 



Bond et al. 



2005 



Dawson et al. 



2006 



scales as cr|'^ (IKomatsu fc Seljak. 



2002 



Sievers et al. 



20091 ). the SZ effect could only explain 



around the 44.5 % of the r nost likely B I MA y alue (14.88 fiK). The Cosmic Background 



Imager (CBI) observations^ 



and, in a recent analysis (jSievers et al. 



Bond et al. 



(|2005r) . cover the multipole range 400 < i < 4000 



20091 ). it is argued that, if the CBI excess power 



observed at very small angular scales is explained as a SZ effect for a certain ag value, then 
a power close to 10 fiK, marginally compatible with BIMA measurements at i = 5237, may 
be explained with the same value of ag- 



The level of power we predict for the weak lensing {AWL plus BWL components) 
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may thus be a missing link, simultaneously ensuring the compatibility of WMAP, CBI and 
BIMA observations. We have predicted a lensing effect of a few micro-Kelvin for these 
large angular scales, and hence a smaller SZ component (a smaller value of ag) could lead 
to simultaneous compatibility with WMAP, CBI, and BIMA observations. Aside from the 
large errors bars on the BIMA measurement, we are also cautious about making firm claims 
since lensing and SZ effects are essentially produced by the same structure distributions 
(galaxy clusters and sub-structures involving dark matter and baryons) and, consequently, 
these effects must be strongly correlated. This implies that the spectra of these two effects 
must be superposed in an unknown way (not merely added). In practice, this superposition 
might be analyzed in detail with ray-tracing through hydrodynamical simulations (including 
both baryons and dark matter). More work along these lines is clearly necessary. 

6. Conclusions 

Our first conclusion, in agreement with other work, is that PM simulations are 
inefficient for calculating CMB lensing due to strongly nonlinear structures (see Fig. [3] and 
comments in Section \^72\i . While it is not impossible to simulate these effects with a PM 
code the required resolution is so large that the additional short-scale resolution provided 
by AP3M (or other high resolution N-body technique) is far more efficient at capturing the 
lensing effect. Hence our development of a combined parallel AP3M-ray-tracing code is a 
necessary step to estimate the lensing signal at high i. 

Only the AWL effect (see Section |3]) produced by scales smaller than 42h~^ Mpc 
between redshifts Zj„ = 6 and Zend = has been estimated with AP3M simulations. This 
choice is appropriate for the following reasons: (1) all the strongly nonlinear scales are taken 
into account; (2) omitting scales greater than 42h~^ Mpc makes our ray-tracing procedure 
efficient (see Section [3]); (3) the effect produced by scales greater than A2h~^ Mpc {BWL in 
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Section [3]) can be studied by using the linear version of CMBFAST; and, (4) the lensing due 
to scales smaller than 42h~^ Mpc, at redshift z > 6 {CWL in Section [3]), can be computed 
using standard semi-analytical methods implemented in CMBFAST. 

Given that our study has strongly focused on the AWL signal, we exhaustively 
investigated the numerical issues involved when estimating this signal. Our main 
conclusions can be summarized as follows: 



the lensing contribution between z = 0.2 and z = is negligible 

for each RLS, a few realizations suffice to get a good average Ci{LU) spectrum and, 
moreover, each single simulation gives a very good spectrum rather similar to the 
average one 

RLSs, which are essentially identical except for the preferred directions, give similar 
spectra (the ray-tracing procedure has little variance) 

simulation boxes of 102Ah~^ Mpc are not fully necessary for i > 1000; however, for 
i < 1000, these large sizes lead to the most reliable spectra and, moreover, such sizes 
should lead to very good spectra in the £- interval (1000-2000) 

simulations in boxes of 512h~^ Mpc lead to good Ce{LU) spectra for 1000 < i < 7000 

for 2000 < i < 7000, all the RLSs lie in a region of width ~0.5 fxK, indicating that 
the RLSs give consistent estimates of the signal in this range 



the signal in the range 400 < £ < 700C 



that found elsewhere {e.g. 



Das fc Bode 



is 2.0 ± 0.4 /xK, which is ~ 1.4 /iK higher than 



20081) 



Despite some simulation uncertainties, our code and technique have lead to a robust 
estimate of the lensing effect in the ^-interval (4200,7000), where it clearly dominates 
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the primary anisotropy. Moreover, the estimated power is larger than that obtained with 
semi-analytical methods and the CMBFAST code. We thus suggest that the resulting 
value of a few micro-Kelvin may explain the excess power at high £ in the BIMA and CBI 
observation s. This conclusion is supported by recent studies based on the Millennium 



simulation (ICarbone et al. 



20081 ). where the authors hav e reported a small co ntribution 



from nonlinearity at £ ~ 4100. However, the methods of .Carbone et al. 



mm have been 



designed to build all-sky lensed maps, and do not have the resolution necessary to perform 
an accurate estimate of the weak lensing by strongly nonlinear structures in the ^-interval 
where we have found our main effect. 

Our direct estimation of the potential gradients at photon positions using PP corrections 
from the local dark matter particles appears to be the main origin of the difference between 
our results and other research relying either on planes or grid interpolations. However, 
we emphasize that differences only occur at the large £ > 2000 values we have been 
investigating. Our method employs extremely fine time resolution, namely that used by 
AP3M simulations and also a very good angular resolution (see above). Due to current 
limitations in our method it is not possible for us to determine the precise role of temporal 
resolution in an accurate lensing calculation. 

It is well known that on small scales baryons do not follow the dark matter distribution. 
Thus, while we have attempted to be as accurate as possible in this dark matter simulation, 
we are now probing scales where contributions from bary ons are beg i nning to become 



Jing et al 



(120061 ) considered 



significant. An investigation of the impact of baryons by 
two types of simulations, the first with non-radiative baryons while the second included 
dissipation and star formation. They found that for 1000 < / < 10, 000 an effect of between 
1% to 10% on the weak-lensing shear angular power spectrum, as calculated from dark 
matter alone, was possible. The largest difference was produced in the run with dissipation 
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and star formation, where a 10% increase in the (shear) Qs at ^ = 10000 was observed. 
Future work is definitely necessary to determine the impact of this physics on lensing 
statistics. We plan to conduct simulations with baryons and feedback processes both to 
identify its impact on the signal we find, and also to systematically evaluate the combined 
impact of the SZ effect and weak lensing. 

For calculations that are accurate to high I it seems to be necessary to move CMB 
photons through the simulation box while structures are evolving, which ensures the spatial 
gradients arc accurately calculated on the photon positions (test particles). Once the 
N-body code has been modified to compute spatial gradients at particle test positions, 
there are no theoretical or technical reasons to take a temporal resolution different from 
that defined by the simulation time step. From the theoretical point of view, the time step 
resolution is obviously most compatible with the N-body technique (and thus takes into 
account the entire evolution of the simulation) . From the technical point of view, the use of 
the time step resolution requires the minimum memory cost as the number of test particles 
between two successive times, although large, is minimized. Because of these considerations, 
temporal resolution is not a parameter to be varied in our calculations. 

Our AP3M code adapted to CMB lensing calculations can be run for different values 
of the parameters defining the LSs; hence, this code allows us to see how the resulting 
angular power spectra depend on the parameters defining both the N-body simulation and 
the ray-tracing procedure. 
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